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ABSTRACT 

In  this  report  the  relevant,  practical  criteria  required  for  the 
numerical  calculation  of  different  types  of  gas  dynamic  flow  fields  are 
discussed.  A  classification  for  finite  difference  schemes  is  given.   Numerical 
experiments  are  conducted  on  the  non-linear  Burger's  equation  using  the 
Brailovskaya,  Cheng-Allen,  Dufort-Frankel,  Crank- Nicolson  and  Lax-Wendnoff 
finite  difference  schemes  as  typical  examples  of  different  classes  of  finite 
differences  schemes.   The  effect  of  time  and  space  increments,  types  of 
boundary  conditions,  perturbations  at  the  boundaries,  truncation  errors,  con- 
vergence criteria  and  convergence  rates  are  studied.   It  is  shown  that  both 
the  Cheng-Allen  and  Dufort-Frankel  schemes  have  individual  merits  in  different 
types  of  application. 
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I .   INTRODUCTION 

Multidimensional  flow  problems  are  governed  by  physical  laws  which 
are  mathematically  represented  by  a  set  of  partial  differential  equations.   In 
attempting  the  solution  of  a  multidimensional  flow  problem,  one  has  recourse 
to  (a)  experimental  methods,  (b)  pure  analytical  methods,  and  (c)  computational 
methods  using  high  speed  electronic  computers.   Experimental  techniques  are 
usually  confined  to  particular  problems  and  yield  a  limited  amount  of  informa- 
tion, while  pure  analytical  methods  can  be  applied  to  the  most  elementary  cases, 
with  many  simplifying  assumptions.  While  each  of  these  two  methods  have 
their  own  individual  and  peculiar  merits,  computational  methods  with  high 
speed  electronic  digital  computers  of  the  present  and  future  generations, 
like  the  ILLIAC  IV  parallel  processing  machine,  offer  a  number  of  advantages, 
among  which  are  the  following: 

1.  A  numerical  solution  of  the  original  partial  differential 
equations  governing  the  problem  furnishes  insight  and 
valuable  information  about  the  inherent  nature  of  the 
problem  without  any  simplifying  assumptions,  or  with  a 
minimum  of  simplification. 

2.  Once  a  computational  program  has  been  developed,  a  wide 
spectrum  of  results  and  any  required  parametric  studies 
are  obtained  with  the  least  amount  of  effort  and  expendi- 
ture . 

3.  A  single  class  of  numerical  methods  can  often  be  applied 
to  a  number  of  physical  problems  with  different  boundary 
conditions,  whereas  a  physical  experiment  would  necessitate 
expensive  changes  in  equipment. 

h.      Computational  procedures  yield  cross  correlations  between 
physical  variables  that  are  not  readily  available  from 
experimental  techniques. 

5-   Optimization  of  the  related  variables  under  given  con- 
straints and  under  physically  realistic  conditions  is 
best  accomplished  computationally. 

6.   Numerical  techniques  are  the  only  experimental  tools 
available  in  certain  problems  because  the  nature  of 
the  physical  phenomenon  and  its  gigantic  scale  precludes 
any  form  of  physical  experimentation. 
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7.   Coupled  with  experimental  investigations  and  asymptotic  analysis, 
computational  techniques  are  a  very  powerful  tool  in  obtaining 
a  complete  understanding  of  the  physical  problem. 

Physical  phenomena,  in  general,  and  multidimensional  flow  problems,  in 
particular,  can  be  described  mathematically  by  a  set  of  conservation  laws  in 
integral  form.   The  equations  are  appropriately  called  conservation  laws  be- 
cause they  express  the  principle  of  conservation  of  mass,  momentum  and  energy. 
These  integral  relations  are  readily  reducible  to  a  system  of  partial  differ- 
ential equations  governing  the  flow.   Together  with  prescribed  initial  and/or 
boundary  conditions,  the  multidimensional  flow  problem  reduces  to  a  boundary 
value  problem  in  the  steady  case  or  a  coupled  initial -boundary  value  problem' 
in  the  non-steady  case.   Computationally  it  is  oftentimes  easier  to  solve  even 
a  steady  state  problem  as  the  long  time  limit  of  the  corresponding  non-steady 
problem,  since  implicit  difference  schemes  employing  iterative  techniques  are 
generally  required  for  the  steady  state  forumulation  [l].   Moreover,  in  multi- 
dimensional flow  problems,  a  time  dependent  formulation  has  decided  advantages, 
even  if  for  no  other  reason  than  that  it  is  computationally  simple.   Such  a 
method  of  computing  the  steady  state  solution  as  the  long  time  limit  of  a 
sequence  of  temporal  states,  is  often  referred  to  as  the  "asymptotic  method." 

In  the  numerical  solution  of  the  partial  differential  equations  of 
multidimensional  fluid  dynamics,  the  computational  space  is  divided  into 
smaller  meshes  or  cells.   The  partial  derivatives  in  the  governing  equations 
are  approximated  by  finite  differences  of  values  at  the  mesh  points.   The 
problem  is  thus  discretized  and  the  dependent  variables  are  calculated  at  a 
discrete  set  of  mesh  points  comprising  the  computational  grid  in  the  flow 
field.   The  usual  assumption  made  is  that  as  the  computational  grid  becomes 
smaller  and  smaller,  the  computed  solution  at  a  point  converges  to  the  actual 
true  solution. 

Numerical  methods  may  be  divided  into  three  categories,  depending 
on  the  construction  of  the  computational  grid  and  the  calculation  of  the 
dependent  variables.   A  brief  description  of  the  three  methods  serves  to 
illustrate  the  advantages  and  disadvantages  of  each. 
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1.1   The  Method  of  Finite  Differences 

This  widely  used  method  is  applicable  to  both  linear  and  non-linear 
partial  differential  equations  of  the  hyperbolic,  parabolic  and  elliptic 
type.   It  stems  from  the  original  work  of  Courant,  Friedrichs  and  Levy  [2]  • 
The  computational  region  is  subdivided  into  smaller  meshes  which  can  have 
rectangular,  cylindrical,  spherical  or  other  frames  of  reference,  depending 
on  the  nature  of  the  problem  and  its  boundary  (Figure  l).   The  grid  lines 
are  arbitrary. 


The  required  dependent  variables  at  any  grid  point,  fixed  in  time 
and  space,  are  obtained  as  a  function  of  the  variables  at  neighboring  grid 
points.   Computational  ease  and  programming  simplicity  are  greatly  facilitated 
if  no  moving  or  stationary  discontinuities  exist  in  the  flow.   Whenever  dis- 
continuities do  exist  in  the  flow,  special  procedures  are  required  to  handle 
them  depending  on  whether  they  are  moving  discontinuities  or  stationary  dis- 
continuities . 

(a)  Stationary  Discontinuities 

When  the  number  of  singularities  is  small  and  the  discontinuities 
are  fixed  in  time  and  space,  the  computational  space  can  be  treated  as  being 
composed  of  a  number  of  separate  regions .   These  regions  have  continuous 
solutions  bounded  by  the  singular  discontinuous  surfaces  which  satisfy  boundary 
conditions  appropriate  to  the  problem.   Finite  difference  methods  are  used  in 
the  regions  with  continuous  solutions . 

(b )  Moving  Discontinuities 

(l)   The  initial  formation  and  subsequent  motion  of  the  discontinuity 
may  be  followed  by  keeping  track  of  the  characteristics  of  the  flow,  which 
are  the  lines  of  propagation  of  discontinuities.   Once  the  discontinuity  has 
been  thus  located  in  space  and  time,  it  may  b©  treated  as  a  boundary  point 
or  surface,  and  the  method  of  finite  differences  used  for  calculating  the 
values  at  the  remaining  points . 
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FIGURE    1(a)     ARBITRARY    GRID    IN    THE    METHOD  OF    FINITE    DIFFERENCES 


CHARACTERISTIC 
LINES 


FIGURE    Kb)     GRID   FORMED  BY   CHARACTERISTIC   LINES 
IN    METHOD    OF   CHARACTERISTICS 


FIGURE    1(c)    CURVILINEAR    GRID   CONFORMING    TO   BODY   SHAPE    USED 
IN    METHOD    OF    INTEGRAL    RELATIONS 
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(2)  In  order  to  avoid  the  additional  complication  in  programming  logic 
as  a  result  of  finding  the  location  of  the  moving  discontinuity,  certain 
devices  may  be  employed  which  permit  the  use  of  finite  differences  alone  through- 
out the  flow,  without  the  necessity  of  resorting  to  other  methods  like  the 
method  of  characteristics.   These  devices  are  based  on  the  physical  nature  of 
the  discontinuous  phenomenon,  as  for  example,  the  use  of  artificial  viscosity 
for  smearing  out  shocks  in  hydrodynamic  flow  [3].  Such  devices  for  computing 
across  discontinuities  should  be  used  with  caution  as  they  introduce  errors 
into  the  computation,  resulting  in  loss  of  accuracy  and  resolution. 

(3)  Where  feasible,  problems  with  non-stationary  singularities  may  be 
solved  by  using  a  movable  grid  attached  to  the  singularity.   This  involves 
prior  knowledge  of  the  motion  of  the  discontinuity.   The  transformation 
should,  however,  avoid  the  introduction  of  independent  variables  which  map 
such  grids  into  coordinate  lines. 

1.2  The  Method  of  Characteristics 


Numerical  solutions  of  partial  differential  equations  by  the  method 
of  characteristics  are  possible  only  for  equations  of  the  hyperbolic  type. 
In  the  method  of  finite  differences,  the  computational  grid  may  be  constructed 
for  any  reference  system  of  coordinates  and  may  be  superposed  on  the  com- 
putational space  at  the  very  outset.   However,  in  the  method  of  characteristics, 
the  computational  grid  must  be  bounded  only  by  characteristic  lines  or  suf- 
faces  which  are  drawn  as  the  calculation  proceeds  in  space  and  time.   (See 
Figure  lb.)   This  is  necessary  since  it  is  only  along  these  characteristic 
lines  or  surfaces  that  the  original  non-linear  partial  differential  equations 
can  be  reduced  to  simpler  ordinary  differential  equations,  called  the  com- 
patibility conditions.   In  the  three  dimensional  case,  complexity  of  computa- 
tional algorithm  and  difficulties  in  programming  arise  because  of  the  complex 
behavior  of  the  characteristic  surfaces. 

The  method  of  characteristics  is  eminently  suitable  for  flows  con- 
taining a  limited  number  of  discontinuities,  as  it  permits  the  accurate  deter- 
mination of  the  origin  and  subsequent  propagation  of  shocks.   Also,  the  method 
admits  of  considerable  mathematical  vigor,  and  uniqueness  and  convergence  have 
been  proved. 
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1. 3  The  Method  of  Integral  Relations 

This  method,  like  the  finite  difference  method,  is  applicable  to 
hyperbolic,  parabolic  and  elliptic  partial  differential  equations.   The  com- 
putational space  is  divided  into  a  series  of  curvilinear  strips  conforming  to 
the  shape  of  the  boundary  of  the  region  of  interest.   (See  Figure  lc.  )   The 
governing  partial  differential  equations  are  written  in  divergence  or  con- 
servation form  and  are  integrated  .across  these  strips,  with  the  functions  occur- 
ring in  the  integrands  being  replaced  by  some  interpolation  expressions. 
Thus,  in  effect,  it  is  the  integrals  that  are  being  approximated,  instead  of 
the  partial  derivatives,  as  in  the  finite  difference  method.   Since  the 
governing  equations  are  integrated  in  divergence  or  conservation  form,  the 
results  obtained  therefrom  remain  valid  even  when  crossing  a  discontinuity. 

Each  of  the  above  methods  may  be  applied  individually  or  in  conjunc- 
tion with  any  of  the  others  depending  on  the  nature  and  complexity  of  the 
problem.   To  name  two  examples,  the  supersonic  flow  over  a  satellite  nose  cone 
containing  an  oblative  shield  may  be  calculated  by  a  combination  of  the  method 
of  characteristics  and  the  method  of  integral  relations.  Flow  of  a  viscous 
fluid  in  a  boundary  layer,  may  be  calculated  by  the  method  of  finite  differ- 
ences, the  method  of  integral  relations,  or  both. 

Another  very  important  factor  to  be  considered  in  the  selection 
of  any  one  of  the  above  three  numerical  methods  is  the  type  of  computer  employed 
and  the  programming  language  used.   For  computers  like  the  parallel  processing 
ILLIAC  IV,  maximum  advantage  must  be  taken  of  its  capability  to  compute  simul- 
taneously in  parallel.   If  the  method  of  characteristics  were  used,  the  com- 
putational grid  cannot  be  initially  superposed  on  the  flow  field,  and  hence, 
until  the  flow  has  developed  sufficiently,  all  the  parallel  processors  cannot 
be  put  into  operation  simultaneously. 

In  the  method  of  finite  differences,  the  division  of  the  computa- 
tional space  into  arbitrary  grids  can  be  accomplished  at  the  very  beginning. 
In  many  problems,  calculations  at  several  grid  points  proceed  simulta- 
neously and  independently,  so  that  each  processor  or  arithmetic  unit  of  a 
parallel  processing  machine  can  be  assigned  to  a  given  grid  point.   In  this 
way  great  saving  in  computational  time  is  realized,  making  possible  the  cal- 
culation of  solutions  to  problems  that  would  otherwise  be  prohibitively  time 
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consuming.   In  the  parallel  processing  machine,  the  computation  of  points  ad- 
jacent to  shock  discontinuities  and  of  boundary  points  presents  certain  chal- 
lenges.  In  the  computation  of  shock  boundaries,  one  simple  expedient  would  be 
to  incorporate  a  system  of  checks  which  would  he  processed  by  all  the  arith- 
metic units,  but  would  be  effective  only  if  shocks  were  present. 

Since  the  method  of  finite  differences  offers  definite  advantages, 
especially  from  the  viewpoint  of  application  to  advanced  parallel  processing 
machines,  we  shall  confine  our  attention  to  this  numerical  technique.   For  com- 
putation of  three  dimensional  flow  problems,  finite  difference  techniques  offer 
the  maximum  assistance  in  the  natural  extension  from  the  two  dimensional  case. 


-T- 


2.   CONSERVATION  LAWS 

The  laws  that  govern  the  motion  of  a  fluid  or  gas  state  that  given 
a  control  volume  V  bounded  by  a  surface  S,  the  sum  total  of  all  the  phenomena 
occurring  within  the  volume  V  must  be  such  that  the  overall  mass,  momentum  and 
energy  contained  within  the  control  volume  are  conserved.   Expressed  mathemat- 
ically, the  conservation  laws  may.be  written  as  follows: 

Conservation  of  mass 


where 


_9_ 

at 


«/v  pdV  =  /s   "   P^  '  ^ 


(1) 


Conservation  of  momentum 


3t 


/   pu  dv  =  9  -  (pu)  u  -  pi  +  t   .   ds  +  I  pg  dv  (2) 


Conservation  of  energy 


9t       J  v 


L 


u 
2 


p(e  +  —  )   dv 


=    fs  ^   pu 


+  \    pHdv+L 


(e  +  £)   +  u 


T    -   U      •     lp 


] 


.    Ip      .    ds 


pHdv  +  J      pg    .   u  dv 


(3) 


u  =  velocity  vector 

p  =  density 

p  =  pressure 

-*■ 

t  =  shear  stress  tensor 


->■ 
-»■ 
I  =  unit  tensor 

H  =  heat  conduction  from  external  sources 

g  =  body  force  vector 

While  a  number  of  solutions  will  satisfy  the  above  conservation  laws, 
the  particular  solution  for  any  given  problem  is  determined  uniquely  by  its 
initial  and  boundary  conditions.   The  choice  of  a  reference  coordinate  system 
in  a  computational  procedure  employing  finite  differences  is  governed  by  the 
ease  of  obtaining  a  finite  difference  formulation  as  well  as  suitability  for 
programming  in  various  languages. 

In  particular,  in  Cartesian  tensor  notation,  the  conservation  equa- 
tions may  be  written  as  follows: 


Conservation  of  mass 


&  ♦ -i- (pu  )  -  o  w 


J 


Conservation  of  momentum 


|t-  (pu,)+r^-  (pu.u  -)+!£-«.,  -■—-  x.  -pg.=0  (5) 

dt      1     dX.       1  J     9x.   1J    3x.    1J       1 

J         J      J 


Conservation  of  energy 


ie. 
at 


2  3     i-  2  -j 

^   +  ^    +  Si     "    ^   +  ^    +   PUJ    "  Ui    'ij  J 


-   Pgi  u±   +      pH  =  0  (6 


where   &..    is   the   Kronecker   delta. 
1J 
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The  equations  (k)   through  (6)  can  readily  he  obtained  from  the  in- 
tegral equations  (l)  through  (3)  "by  applying  Gauss  theorem 

Js   b.  n  ds  =  Jv     V.  b  dV  (7) 

2.1  The  Navier-Stokes  Equations   - 

If  we  assume  that  Stokes  law  of  fluid  friction  applies,  i.e.,  the 
shear  stress  is  proportional  to  the  rate  of  shear  deformation,  and  if  we  assume 
that  the  constant  of  proportionality  or  viscosity  coefficient  y  is  a  constant, 
and  that  the  bulk  viscosity  coefficient  is  zero,  we  may  write  the  shear  stresses 
as  follows  for  the  two-dimensional  case: 


t^  -  ~    y  —  -  o  v  —  (8) 


k 

=  3     y 

3u 
3x  " 

2        3v 

3VW 

h 
=  3  y 

3v 

ay 

2  3u 

3  y  ^ 

v"'w   ---"-  (9> 


/3u   3v\  /,«\ 

xy    yx       3y   3x 


When  these  values  for  the  shear  stress  t. ,  are  incorporated  into  the  conser- 
vation  equations,  they  are  called  the  Navier-Stokes  equations. 

It  will  now  be  shown  that  the  conservation  laws  embodying  the 
Navier-Stokes  shear  stress  has  the  same  vector  form  as  the  Burgers '  equa- 
tion [h]: 


3u     3u  =  1_  afu  (11) 

3t   U  3x   R  .  2 
e  3x 
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In  the  absence  of  the  external  body  force  g  and  heat  addition  from 
external  sources,  (which  does  not  include  heat  generated  or  absorbed  within 
the  body  from  such  causes  as  no n- equilibrium  chemical  processes,  etc.)  we  may 
write  the  conservation  laws   in  one  dimension,  thus: 


Conservation  of  mass 


£+|jCp»)-0  (12) 


Conservation  of  momentum 


ft  <>»> +  k (pu2  +  p  - 1  *  &  ■  °  (13) 


Conservation  of  energy 


2  2 

3t    (e  +   2    }    +  ix"     [pu(e  +  2)+pU'U'3^]=°  {lk) 


If  we  now  define  the  vectors 


(15) 


-  pu 

■(pu2  +  p) 

i-  pu2  (p+E)j 
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•-'W-2  I  (17) 

ax"(u  3  y^7 


where  E  =  p  (e  +  |-)  (l8) 


we  can  write  [5]  one  vector  equation  for  the  three  conservation  laws  (12), 
(13),  and  (lU)  as 


U=  A  U  +  B  U  (19) 

t        X        XX 


where         A  =  —  and  B  =  — 

x 


Note  that  the  conservation  laws  have  the  same  vector  form  as  the  Burgers' 
equation  (ll) 


2 
9u  .  _  _9u   1_  9  u 

3t  :=  "U  3x   R   .2 
e   9x 


The  Burgers'  equation  (ll),  therefore,  serves  as  a  useful  model  for  testing 
finite  difference  schemes,  prior  to  actual  use  on  the  full  Navier-Stokes  equa- 
tions.  In  this  way,  flexibility,  ease  of  programming  and  a  substantial  saving 
in  computational  time  are  secured  when  testing  a  number  of  difference  schemes, 
with  the  aim  of  studying  the  different  characteristics  of  each  scheme. 
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3.   DIFFERENCE  SCHEMES 


The  numerical  solution  of  a  given  set  of  partial  differential  equa- 
tions involves  replacing  the  partial  derivatives  by  finite  differences  which 
are  functions  of  the  discrete  values  at  neighboring  points.   The  finite 
difference  representation  is  often  involved  and  bears  no  resemblance  to  the 
original  partial  derivative.   The  choice  of  a  difference  scheme  for  a  given 
problem  involves  careful  consideration  of  a  number  of  factors*  including 
the  following: 

1.  Consistency  of  the  difference  approximation. 

2.  Convergence  and  stability  of  computed  results. 

3.  Truncation  errors  and  accuracy. 
h.  Ease  of  programming. 

5.  Reasonable  computation  time. 

6.  Adaptability  to  different  programming  languages. 

7.  Applicability  to  present  and  future  generation  of  computers. 

A  brief  description  of  various  types  of  difference  schemes  and 
various  practical  considerations  pertaining  to  their  use  will  be  given  in 
this  section  and  in  Section  h. 

3.1   Types  of  Difference  Schemes 

The  one-dimensional  Burgers1  equation 


2 
9u  L    3u_  _  _1   9  u 

3t   U  3x  "  R   .  2 
e   3x 


may  also  be  represented  as 


u  =  -  f(u)  u  +  |-  u  (20) 

t  x   R    xx 

e 


where  f(u)  =  u 
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Let  u  be  the  value  of  the  computed  solution  at  point  (j,n)  in  com- 


putational space  (Figure  2). 


3-2   j-1   j     j+1   j+2 


n+2 

n+1 
n 

n-1 
n-2 


Figure  2 


The  values  at  neighboring  points  u..-,   ^>«  and  u  .—  '  ..  can  be  obtained  as 

j+1,  2        j 

functions  of  the  value  u  .   at  point  (j,n)  thus, 

J 


u  .   =  u  .  +  At  u,      +  0  ( At ) 


(21) 


In  this  representation,  we  have  made  the  following  assumptions: 

(a)  The  Taylor  expansion  about  the  grid  point  (j,n) 
converges  mo not on ic ally . 

(b)  This  convergence  is  rapid  enough  to  permit  truncation 
of  the  series  after  the  second  term. 


Therefore,  9ui 
3t  I 


.n+l   n 
u .    u  . 

J  aI  ''   +  °  (it) 


Thus    forward   (and  backward)     differencing  involves   an  error  of  first   order   in 
(At). 
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The  second  term  in  the  Burgers  equation  (11 )  can  be  represented  as 


and  the  third  term  as 


1     _  1 r  m        p    q   i 

5—  u   =  Lu.   -  2  u .  +  u .  _  J 


e        Ax^R 


e 


where  k,   Z,   m,    p  and  q  are   either  n+1   or  n.      Thus   the  Burgers    equation  may  "be 
represented  in  general   form  as: 


n+1   _     n       At        n   ,   k               I      ,        At        r   m  _     p  q      , 

u^        =  u^   "  ^TZ  f<    vu.^    -  u,   n  )   +  [u,^    -  2  <  +  u*  ±] 

(22) 


j  j  2Ax      j  j+1  J-l  A    2  J+1  J  J-l 


(a)      Explicit   Schemes   -  One  step  schemes 

If  k  =   £  =  m  =  p  =   q=n,   we  have  the  Euler's   scheme 


n+1  n       At      ji    /    n  n      v        At        /    n  „  n  n      x 

Uj        =  Uj    "  2^  fj    (Vl   "  Uj-l'    +  ^  (Vl   "  2UJ    +  Uj-l' 

(23) 


This  is  an  explicit  scheme  where  values  at  the  (n+l)   time  step  are  calculated 
directly  from  known  values  at  the  n   time  step.   No  intermediate  or  provisional 
values  of  u    are  calculated. 

J 

(b)      Explicit  Schemes   -  Two   step  schemes 

Quite  often  two-step  explicit   schemes   are  encountered  of  which   the 
Richtmyer  variation  of  the  Lax-Wendroff  scheme   is    a  familiar  example    [6], 
In  the   inviscid  case    (R     ■>  °°),  the  scheme  may  be  written  as 
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in  n 

un+2   =  _4+l I  +  At_f       Itu  -   u    ] 


1  1  1_ 

2  r    n+2  n+2 


n+1  n       At       .n+2    r    n+2  n+2,  /„i  x 

u.        =u.   +^     f  [u         -  u  j]  (2k) 


In  this   scheme,    intermediate   values  centered  at  time   (n+y)  are   first  obtained 
before  the   final   values   at  time    (n+l)    are   calculated.      The  intermediate  values 
may  also  be   calculated  at  time    (n+l)    instead  of   (n+y)  as    is    done   in  the 
Brailovskaya      and  Cheng-Allen  schemes   to  be   discussed  later. 

(c)  Fully  Implicit  Schemes 

If  k  =    £  =  m  =  p  =  q  =    ..    n+l,    the  calculation  of  u.   "   is   more 

J 

involved  and  iteration  techniques  are  necessary.   The  fully  implicit  schemes 
can  be  shown  to  be  unconditionally  stable. 

(d)  Partially  Implicit  Schemes 

If  either  the  convective  term  or  the  viscous  term  alone  is  formu- 
lated implicitly  and  the  other  term  is  formulated  explicitly,  we  have  a 
partially  implicit  scheme.   Thus,  the  convective  term  is  formulated  implicitly 
if 

k  =  I   =  n+l  and  m  =  p  =  q  =  n, 


and  If 


k  =  I  =   n  and  m  =  p  =  q  =  n+l 


n+l 
the  viscous  term  is  formulated  implicitly.   Calculation  of  u.   again  involves 

J 

iteration  techniques. 
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(e)   Quasi  Implicit  Schemes 

If  k  =  m  =  n  and  £  =  p  =  q  =  n+1, 


n+1    n   At   ji/  n      n+1,.    At 
u.   =  u.  -  ttt—  f  .(u.  ...  -  u.  n  )    +  — r=- 
J      j    2  Ax  j   j +1    j-1    Ax2R^ 


n+1 


r  n     _  n+1    n+1 -■ 
u.    -  2u.   +  u,   J 

(25) 


Values  u.  "  can  be  calculated  explicitly  without  any  iteration 

J   n+1  . 
techniques  provided  u    is  known  at  the  left  boundary.   The  scheme  is  then 

called  a  quasi-implicit  scheme. 


( f )   Differencing  in  Integral  Form 

Finite  difference  schemes  may  be  obtained  directly  from  the  integral 
conservation  laws  expressed  in  Equations  (l)  -  (3).   This  method  has  been  used 
by  Roberts  and  Weiss  [7]  for  unsteady  convective  problems  and  possesses  certain 
advantages,  especially  in  flows  containing  discontinuities,  since  the  form  of 
differencing  embodies  the  conservation  laws. 

In  a  two-dimensional  problem,  a  computational  net  may  be  defined  as 
in  Figure  3.  An  elemental  volume  AV  is  then  given  by 


AV  =  Ax 


Ay  =  (x+  -  x_)  (y+  -  y_) 


x 


-y. 


Figure  3 


Using  this  difference  formulation,  Equation  (l)  ''may  be  written  as 
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3       f  f  y+  I  X+  f  X+ 

—  JAVpdxdy  =  J  -    (pu)    dx  I  +J  x      -    (pu)dx 


(26) 


and  similarly  for  Equations  (2)  and  (3). 

Another  classification  used  in  addition  to  the  above  is  based  on 
the  number  of  time  levels  used  in  the  difference  formulation.   Thus,  the 
Dufort-Frankel  and  leap  frog  scheme  can  be  calculated  like  an  explicit  scheme, 
but  uses  data  points  from  two  time  levels  n  and  (n-l)  to  calculate  points  on 
the  time  level  (n+l).   It  is,  therefore,  called  a  "two  level"  scheme.   Multi- 
level schemes  vary  from  the  fully  explicit  to  the  fully  implicit  category. 
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CONSISTENCY,  ERRORS,  STABILITY  AND  CONVERGENCE 
OF  DIFFERENCE  SCHEMES 


k .1   Consistency 

In  finite  difference  methods,  the  continuous  partial  derivatives 
are  approximated  by  functions  or  operators  involving  known  values  at  a 
selected  number  of  discrete  points.   The  solution  marches  on  from  a  given 
sequence  of  known  initial  and/or  "boundary  values  to  another  sequence  of  com- 
puted values,  and  so  on,  until  the  final  solution  is  realized,  also  at  a  set 
of  discrete  points. 

The  difference  formulation  is  said  to  be  a  consistent  approximation 
if  the  final  computed  solution  represents  the  true  solution  to  the  problem 
at  the  selected  discrete  points. 

To  insure  consistency,  a  number  of  important  practical  guidelines 
must  be  observed. 

(l)   Since  the  finite  difference  computation  is  a  step-by-step 
marching  procedure,  the  difference  scheme  must  be  consistent  at  each  step 
of  the  calculation. 

Consider  the  calculation  of  u.   from  the  known  value  u1?  for 
-k  <_  j  <_  k,  where  k  is  a  positive  integer.   Further  assume  that  U.   are 

J 

exact  values.   If  0  (At)  is  a  finite  difference  operator  such  that  it  trans- 

_       .  .    n  ,    n+1 
forms  points  u.  to  u.   ,  i.e., 
J     J 


n+1 


0  (At)  un  =  u 
J    J 


then  the  difference  operator  0  (At)  provides  a  consistent  approximation  if 


Un+1  -  un+1  |  |  +   0  as  At  -*  0,  0  <  t  <_  T  (27) 

J  J 


Consistency  of  the  overall  calculation  is  secured  by  insuring  consistency  at 
each  step  of  the  calculation. 
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(2)   The  quantity  under  the  norm  in  (27)  can  be  written  as 


U  (t  +  At)  -  0  (At)  u(t; 


This,  in  effect,  is  equivalent  to  saying  that  the  truncation  error  of  the 
difference  operator  0  (At)  must  tend  to  zero  as   At  -*•  0. 

(3)   Since  the  solution  marches  from  known  values  to  a  set  of  com- 
puted values,  it  would  initially  begin  with  a  known  sequence  of  exact  bound- 
ary values.   To  insure  that  the  final  computed  solution  is  consistent  with 
the  true  final  solution,  it  is  important  that  the  finite  difference  formula- 
tion of  the  boundary  conditions  be  consistent  with  the  true  physical  condi- 
tions at  the  boundary. 

(h)      The  final  computed  solution  should  depend  continuously  on  the 
given  initial  and/or  boundary  conditions,  once  these  boundary  conditions  are 
consistently  represented.   Since  disturbances,  in  the  form  of  changes  in  the 
flow  field,  propagate  along  characteristic  lines  ,  the  final  solution  will 
depend  continuously  on  the  initial  data  if  all  points  in  the  computational 
space  are  calculated  from  known  data  within  the  domain  of  dependence  of  the 
point  being  calculated.   This  principle  is  embodied  in  the  Courant-Friedrichs- 
Lewy  rule  which  states  that  the  time  increment  At  should  satisfy  the  relation 


At  <    ,AX  (28) 


u  +  a 


where  u  is  the  flow  velocity  and  a  the  speed  of  sound  in  the  gas . 

This  condition  also  states  that  discontinuities  in  the  flow  should 
be  properly  demarcated,  so  that  the  characteristic  lines  do  not  cross  them. 
Otherwise,  the  final  solution  would  not  depend  continuously  on  the  initial  data, 
This  precludes  the  possibility  that  the  operator  0  (At)  would  contain  points 
on  the  opposite  side  of  a  discontinuity. 


-20- 


h .2     Errors    and  Accuracy 

Consider  the  equation 


^  +  u^=0  (29) 

3t  3x 


and  a  finite  difference  formulation  to  it 


n+1    n      /  n      riv 

U.     -  U.         (U.  ...  -  U.  J 


L(u)  -  '1   At   J  ♦  u^   '1+1Ax   J   =  0  (30) 


Expanding  u.   and  u    about  the  point  u . f   using  Taylor's  series  in  powers  of 
t  and  x  respectively,  it  is  shown  that 


2  2 

L(U)  "  (9?+  U^}  =T  9t2  +  ••'  +—  9xT+  •••        (31) 


The  error  of  truncating  the  Taylor  series  is,  therefore,  of  order  At  and  Ax. 

Besides  truncation  errors,  other  types  of  errors  introduced  into  a 
calculation  are: 

(a)  Boundary  errors  introduced  into  the  computational 
space  at  the  boundary,  by  incorrect  formulation  of 
the  boundary  conditions. 

(b)  Inconsistency  errors. 

(c)  Machine  roundoff  errors. 

These  errors  propagate  and  accumulate  in  complex  fashion  to  yield  a  gross 
error  in  the  computed  result. 
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k.3     Stability  and  Convergence 

The  finite  difference  solution  to  a  given  system  of  partial  differ- 
ential equations  involves  a  marching  procedure  applied  to  a  sequence  of  dis- 
crete points  in  computational  space.   The  difference  operator  0  (At)  applied 

once  transforms  initial  values  u.  to  values  u.  at  the  first  level  of  com- 

J  J 

putation.   Let  N  levels  of  computation  be  required  to  attain  the  final  solu- 
tion.  Since  the  operator  0  (fat)   contains  an  intrinsic  error,  there  exists 
the  possibility  that  the  cumulative  error  over  N  computational  levels  becomes 

unmanageable.   Hence,  the  computation  is  said  to  be  stable  if  given  a  value 

j 
At  such  that  0   <_  At  <_  t,  0  (At )   is  uniformly  bounded  for  0  <_  N  At  <  T, 

where  T  is  the  time  at  the  final  solution. 

This  final  solution  is  said  to  converge  to  the  true  solution  at  any 
time  level  if  the  norm  of  the  difference  between  the  computed  solution  and  the 
true  solution  tends  to  zero  as  At  tends  to  zero. 

Let  U.  be  the  true  solution  at  point  (j,n),  and  u  be  initial  values. 
J  o 

The  difference  approximation  is  convergent  if 


0  (At)n  u  -  Un   I  ->  0  as  At  ->  0  (32) 

o    J  '  ■ 


for  all  j  . 

Stability  and  convergence  are  synonymous  as  proved  by  P.  Lax  for 
linear  partial  differential  equations.   Lax's  equivalence  theorem  [8]  states 
that,  given  a  properly  posed  initial  value  problem  and  a  consistent  finite 
difference  approximation  to  it,  stability  is  the  necessary  and  sufficient  con- 
dition for  convergence.   Although  the  theorem  is  proved  only  in  the  linear 
case,  it  is  quite  common  to  assume  that  it  applies  locally  to  the  quasi- 
linear  partial  differential  equations  of  fluid  dynamics . 

In  practical  computations,  the  non-linear  partial  differential 
equations  do  not  readily  admit  to  a  simple  stability  analysis.   At  best, 
stability  is  investigated  with  a  simplified  linear  form  of  the  original  equa- 
tions.  The  Von-Neumann  analysis,  which  is  the  most  commonly  used  procedure 
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for  testing  stability,  studies  the  behavior  of  a  sinusoidal  error  superposed 
on  the  true  solution.   If  the  amplitude  of  the  error  increases  as  calculation 
proceeds,  the  difference  operator  is  considered  unstable.   A  more  difficult 
means  of  studying  the  stability  of  a  difference  operator  is  the  so-called 
energy  method,  which  provides  information  at  both  interior  and  boundary  points. 
The  essential  feature  of  this  method  is  to  show  that  the  increase  in  a  specific 
norm  of  the  solution  vector  is  no  greater  than  1  +0  (At)  between  the  two  suc- 
cessive time  steps.   If  this  norm  can  then  be  shown  equivalent  to  the  L2  norm, 
stability  in  this  norm  is  assured.   In  fairly  simple  cases,  the  physical 
energy  of  the  system  provides  such  a  norm,  and  hence  the  method  is  called  the 
energy  method.   The  mathematical  analysis  is  quite  complicated,  and,  unlike 
the  Von- Neumann  analysis,  the  energy  method  does  not  provide  a  physical 
insight  into  the  manner  of  error  propagation. 

In  practice  the  norm  used  in  Equation  (32)  to  test  for  convergence 

takes  various  forms.   Let  u.  be  the  computed  points  and  U.  the  values  of  the 

J  J 

true  solution  at  these  points. 

Let  N  =  a  certain  number  of  iterations 
o 

e  =  a  given  positive  quantity  independent  of  N 
If  the  condition 


TTn    n  i  ,   . 

U.  -  u.    :  e  (33) 


is  satisfied  for  N  >  N  ,  convergence  has  been  attained  in  the  maximum  modulus 
norm.   The  computed  solution  converges  uniformly  to  the  true  solution.   Since 
a  convergence  criterion  based  on  the  maximum  modulus  norm  implies  uniform  con- 
vergence, a  necessary  condition  for  such  a  criterion  to  be  used  is  that  the 

true  solution  U.  must  be  continuous. 
J 

Another  common  norm  employed  to  test  for  convergence  in  numerical 
applications  is  the  mean  square  error  norm,  which  can  be  stated  as  follows: 
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\i     (U?  -  un)2  At   <   e  (3M 

J  J=1     J  J 


for      N>N,0<At<x. 
o  —         — 


This  is  a  measure  of  the  mean  quadratic  error  and  insures  convergence 

in  the  mean.   If  the  sequence  of  points  U.  of  the  true  solution  is  continuous 

J 
and  has  continuous  partial  derivatives,  then  convergence  in  the  mean  and 

uniform  convergence  are  equivalent. 
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5.   NUMERICAL  EXPERIMENTS  WITH  FINITE 
DIFFERENCE  SCHEMES 


In  the  foregoing  the  salient  features  of  obtaining  the  numerical 
solution  of  a  set  of  partial  differential  equations  have  been  described,  and 
among  the  three  main  methods  of  approach,  the  method  of  finite  differences 
has  been  singled  out  for  particular  attention.   The  integral  conservation  laws 
which  are  the  underlying  physical  basis  for  the  governing  partial  differential 
equations  are  given,  and  methods  of  obtaining  their  finite  difference  formula- 
tion are  discussed.   Various  considerations  of  practical  relevance,  in  regard 
to  consistency,  error  propagation,  stability  and  convergence  of  difference 
schemes  have  been  given. 

The  selection  of  a  difference'  scheme  for  a  particular  problem  con- 
stitutes a  first  step  in  the  solution  of  a  given  system  of  partial  differential 
equations.   A  number  of  workers  have  made  comparative  studies  of  different  classes 
of  difference  schemes.   Emery  [9]  has  made  an  evaluation  of  five  difference 
schemes  for  one-  and  two-dimensional,  non-steady  inviscid  flow.  They  are: 

(a)  Lax's  centered  time,  forward  space  difference  scheme 

(b)  Rusanov's  [10]  scheme  which  is  an  improvement  of  Lax's 
scheme,  possesses  the  minimum  artificial  viscosity  at 
each  nodal  point,  and  also  weights  the  importance  of 
neighboring  points 

(c)  Landshoff's  scheme  which  introduces  an  artificial 
pressure  on  the  lines  of  Richtmyer's  and  Von-Neumann's 
artificial  viscosity 

(d)  The  Lax-Wendroff  scheme 

(e)  Richtmyer's  variation  of  the  Lax-Wendroff  scheme. 

Applied  to  a  moving  shock  problem,  most  of  the  schemes  show  a  visible  overshoot 
at  the  shock  front.   It  is  found  that  Lax's  scheme  is  easy  to  program  and  of 
good  resolution.   Rusanov's  method  is  more  versatile  and  also  possesses  good 
resolution.   A  modified  version  of  the  Lax-Wendroff  scheme  is  shown  to  yield 
a  resolution  which  is  spatially  three  times  as  great  as  that  obtained  by 
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Lax's  method  for  the  same  number  of  nodal  points.   Computational  times  are 
minimum  in  Lax's  method  and  Richtmyer 's  variation  of  the  Lax-Wendroff  method, 
and  maximum  in  the  Lax-Wendroff  method.   A  great  deal  of  difficulty  in  pro- 
gramming is  experienced  with  the  Lax-Wendroff  scheme,  even  in  the  inviscid 
case. 

Fourth  order  accurate,  conservative  difference  schemes  derived  directly 
from  the  integral  formulation  of  the  conservation  laws  were  studied  "by 
Roberts  and  Weiss  [7].   The  schemes  are  shown  to  be  free  from  numerical  dis- 
persion, and  thereby  to  yield  accurate  results  in  the  convective  problems 
studied.   The  schemes  are  correctly  centered  in  time  and  space  so  that  the 
modulus  of  the  amplification  factor  is  exactly  unity,  and,  hence,  stability  is 
assured.   Extension  of  the  scheme  to  three  space  dimensions  is  indicated. 

A  similar  technique  based  on  Godunov's  [ll]  method  was  utilized 
by  Taylor  and  Masson  [12]  to  calculate  flow  over  bodies  with  boundaries  of 
large  curvature.   A  curvilinear  co-ordinate  system  was  used.   This  method, 
like  that  of  Roberts  and  Weiss  involves  the  calculation  of  fluxes  across  cell 
boundaries.   The  difference  equation  for  each  cell  is  derived  by  integrating 
the  governing  equations  over  the  volume  of  the  cell. 

Rubin  and  Burstein  [5]  investigated  the  use  of  the  Richtmyer 
two-step  variation  of  the  Lax-Wendroff  scheme  for  both  inviscid  and  viscous 
flow  in  the  case  of  a  moving  shock.   The  full  Navier-Stokes  equations  in  one 
dimension  are  used  in  the  viscous  flow  calculation.   The  difference  formulation 
in  the  viscous  case  is  very  susceptible  to  numerical  instability.   Once  a 
stable  formulation  has  been  successfully  established,  the  Courant-Friedrichs- 
Lewy  stability  criterion  must  be  strictly  enforced. 

More  recently  Rubin  and  Preiser  [13]  generalized  the  Richtmyer 
method  to  three  space  dimensions  and  time.   The  resulting  formulation  is 
obtained  directly  from  the  integral  conservation  laws,  with  some  approximation. 

The  principle  of  the  Lax-Wendroff  scheme  [lU]  was  applied,  in 
non-conservation  form  by  Moretti  and  Salas  [15 ]  to  the  problem  of  shock  forma- 
tion caused  by  an  accelerating  piston.   Consistency  of  formulation  is  maintained 
in  a  transformed  co-ordinate  system,  which  insures  virtually  no  overshoot  in 
the  high  pressure  side  of  the  shock.   Since  the  equations  are  not  maintained 
in  conservation  form,  the  Jacobian  of  the  transformation,  encountered  while 


-26- 


evaluating  the  time  derivatives  in  terms  of  space  derivatives,  is  no  longer 
present.   The  flow  variables  are  evaluated  only  at  the  specified  nodal  points 
and  no  half  incremental  values  are  necessary  as  in  the  Richtmyer  variation 
of  the  Lax-Wendroff  scheme.   Moretti  [l6]  has  also  tested  a  number  of  difference 
schemes  applied  to  the  one- dimensional  equations  in  divergence  form  and  in  non- 
divergence  form. 

A  modification  of  a  scheme  originally  proposed  by  Brailovskaya  [IT] 
was  used  by  Cheng  and  Allen  [l8j .   Stability  of  the  modified  scheme  is 
shown  to  be  independent  of  the  Reynolds  number  for  the  linear  case. 

While  many  different  schemes  were  used  for  specific  problems,  and 
in  some  instances,  comparisons  were  made  among  different  schemes  of  similar 
construction,  a  number  of  questions  still  need  to  be  answered,  among  which  are 
the  following: 

1.  What  is  the  effect  of  the  space  increment  Ax  on  accuracy, 
convergence  rate  and  resolution? 

2.  What  is  the  influence  of  the  time  step  At  on  accuracy  and 
computational  time  in  seeking  a  steady  state  solution? 

3.  What  is  the  effect  of  different  types  of  boundary  conditions 
on  the  different  schemes? 

h.      How  do  different  formulations  of  the  same  boundary  conditions 
affect  the  problem? 

5.  How  do  the  different  convergence  criteria  compare? 

6.  What  is  the  effect  of  different  classes  of  difference 
schemes  on  computational  time? 

7.  What  is  the  effect  of  Reynolds  number  on  accuracy,  stability 
and  convergence  in  a  viscous  flow  situation? 

8.  How  do  the  schemes  tested  compare  with  regard  to  stability, 
accuracy  and  convergence  rates? 

9.  How  do  various  schemes  compare  for  application  to  different 
categories  of  computers,  like  parallel  processing  machines, 
pipe  line  machines,  etc. 

In  order  to  find  answers  to  these  questions,  the  following  finite 
difference  schemes  were  tested  with  the  non-linear  Burgers  equation  [ll]: 
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(l)      Erailovskaya  Scheme 


u^1   _  „» 


At 


u.  .  z 

J   _  -l        r  /,  n      x2        ,    n      ..  -.  1        rn  ^  n         n      ^ 

"     5*t  C(Vl}     "   (UJ-1)]   +  R~A^  [uj+l  "  2uj   +  uj-l] 


_J i  _  zl_  t(u^+n  )      -   (u.    ,  )    ]   +       1        r    n  n  j      n      , 

J_1  R^  [uj+l  "  2uj   +  Uj-1] 

(35) 


2 
Truncation  error  is  of  order  At,  A*  • 


Ax     Ax  Re 


Stability  condition:  At  <  min  (  ^-^ — p^y,  — r —  )  (36) 


max     u 


(2)      Dufort-Frankel  Scheme 


n+1         n-1 
i.        -  u. 

J J 


u.  -    U.  -,  o  o 

-1     r  /  n     \^  /   n     \^n  J-       r    n  n+1 

=  -r~r~  L  (u  •  -,  )     -  u.  J  -i-  — x —  [u.  ,  -  u.       + 

1+Ax  l      j+l;  v    j-l'    J  2       L    j+1  j 


At  °  °   '  Ax  Re 


n  n-1.  ,„\ 

"i-i  "  ui   ]  (37) 


2    2 
Truncation  error:  0(Ax  ,   At  ) 

Ax 

Stability  condition:   At  <  1 — r  (38) 

J  —  max     u  v      ' 


(3)      Cheng-Allen  Scheme 


—n+1         n 

u .         -    u .  -i  0  -1 

<]  J  r  /   n      \2        1    n      \2n  J-        r    n 

=  "  W  [(ui+l}      "    (ui-l)    ]    +  2  [Ui+l 

At  ^  J  J  ReAx2       J   X 


2n?fl    +un  n] 
J  J-1J 
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n+1    n 
-  u . 

a 

At  °  '  ReAx"" 


j      j      1  r  /_n  \2   /_n   \2,     1   r_jn-l   0  n+1   _n+-l. 


(39) 


2 
Truncation  error:   0(At,  Ax  ) 


Ax. 

Stability  condition:      At  < 1 — r  (ko) 

—  max    |u| 

(k)      Crank-Nicolson  Scheme 


n+1  n  /   n      x2        /   n      N2        /   n-KU2        ,   n+lx2 

u.        -  u.  -,  (u.,J      -    (u.    ,)      +-   (u.,,j      -    (u.   n) 

At  ^  2 


n  _n         n  n-i-1        „  n-t-1         nfl 

I  u.n-2u.fu.,+u.,-2u.        +u.n 

f  X  r         J+l  J  J-l  J+l  J Jzi    -I 


Re  Ax 


(in) 


2 
Truncation  error:      0(At,    Ax   ) 

Unconditionally  stable. 


(5)      Lax-Wendroff  Scheme 

Richtmyer's   two   step  method  was  employed  for  this   scheme.      As   in 
the  above  formulations,   bars  denote   the   intermediate   steps. 


_n+l/2         1    /   n  nx        lAt        ,   n      N2        .   n,2    _ 

Uj+l/2  =  -  (UJ.1    +  V    "   2S  t    ^j^      "    <Uj)      ] 


1      At 

+ 


/   n  n  n         n     \ 

(u.,0-u    ,..    -  u .    +-  u.   ,  J 
v    j 4-2  j+1  j  j-l' 
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n      n 


uj  -  ( p }  "  X  Ax"  t  K,l}  "  <Vi)  ] 


At   r   n     ^  n    n   ,        , ,  „  s 
+  ~2—  [  u    -  2u   f  u    ]        (42) 

Ax  Re    J       J    J  ^ 


The  final  scheme  for  calculating  values  at  the  next  time  step  is  given  "by 

n      n     _n+l/2   _n+l/2 

^   .  un  .  *    (un   f  un  >  (  Vl  '  "j-l   Vl/2  "  uj-l/2 

j      j   4Ax   jfl    j-1'  v  2  ; 


A_t   r   n     o.  n    n  ,,  „» 

[  u    -  2u.  4-  u.   ]        (U3) 


Ax2Re    J+1    "J    J"1 


2 
Truncation  error:   0(At,  Ax  ) 


Three  other  formulations  for  the  Lax-Wendroff  scheme  were  tested  and  found 
to  be  unstable. 

According  to  the  classification  in  Section  3,  both  the  Brailovskaya 
and  Cheng-Allen  schemes  are  two-step  explicit  schemes.   The  Dufort-Frankel 
scheme  is  an  explicit,  two- level,  one- step  scheme.   The  Crank-Nicolson  scheme 
is  implicit,  while  the  Lax-Wendroff  scheme  is  a  two-step  explicit 
scheme.   Attention  has  been  focused  mainly  on  explicit  schemes  because  of  the 
ease  in  programming  in  three  dimensions.   The  Crank-Nicolson  implicit  scheme 
is  commonly  employed  in  the  solution  of  the  parabolic  heat  conduction 
equations.   Since  the  Navier-Stokes  equations  are  a  mixed  parabolic -hyperbolic 
system,  we  are  interested  in  observing  the  behavior  of  the  Crank-Nicolson 
scheme  in  comparison  with  the  other  four  schemes. 

Two  sets  of  boundary  conditions  were  chosen  to  test  their  effect  on 
the  final  solution. 


Set  I 


L.  Boundary  u(o)  =  1 
R.  Boundary  u(l)  =  -1 
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Set  II 

L.  Boundary  u(o)  =  1 
R.  Boundary  u(l)  =  0 
for  0  <  x  <  1 

The  exact  solution  of  the  Burgers  equation  for  the  steady  case  is 
given  by 


■Dp 

u(x)  =  A   tanh  [A(B-x)  ^L]  (41+) 


where  A  and  B  are  determined  from  the  boundary  conditions.   For  boundary 

conditions  I,  there  is  a  sharp  discontinuity  at  x  =  0.5>  while  for  B.C. II, 
there  is  a  sharp  discontinuity  at  x  =  1. 

No  special  provision  is  made  for  these  discontinuities  in  the 
numerical  scheme,  to  observe  how  well  they  perform  under  extreme  conditions 
where  discontinuities  may  exist  in  the  flow. 

For  a  meaningful  comparison  of  the  various  schemes,  the  arithmetic 
statements  in  the  computer  program  were  made  as  identical  as  possible,  except 
for  essential  differences  necessitated  by  the  special  characteristics  of  the 
individual  difference  schemes. 

Tests  were  run  at  different  Reynolds  number  Re,  and  at  different 
space  increment  Ax.   The  Reynolds  number  was  varied  from  10  to  1000.   The 
number  of  space  intervals  varied  from  10  to  58  as  60  points  was  near  the 
upper  storage  limit  of  the  CDC  l6ok   computer  used  in  the  calculations.   For 
a  given  space  increment  Ax  and  Reynolds  number  Re,  the  time  step  was  varied 
such  that 

i^<At<Ax  (45) 

k 

This  range  of  At  was  employed  to  subject  the  schemes  to  the  most  difficult 
conditions  possible,  both  at  the  low  Reynolds  number  limit  and  at  the  inviscid 
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high  Reynolds  number  limit.   The  criterion  used  to  test  for  convergence  to 
the  steady  state  was 


max  |unfl  -  un|  <  10" 5  (2+6) 

^    J      J 
J 

at  any  time  step.   This  is  analogous  to  the  maximum  modulus  norm  and  implies 
uniform  convergence.   The  quantity 


1   v   |  /  n+lx2    ,  n>.2i  ,,  . 

—  E   |(u   )   -  (u  )  I  (1+7) 

J  j=l    J        J 


was  also  calculated  for  comparison.   This  is  analogous  to  the  mean  square 

error  norm,  because  u.  "  tends  to  the  true  solution  IT.  as  n  increases.   It 

J  J 

implies  convergence  in  the  mean.   It  is  thus  possible  to  compare  the  two 

norms  for  convergence. 

The  computational  program  is  coupled  to  a  graph-plotting  subroutine, 
which  traces  out  the  required  curves  on  an  oscilloscope  from  which  photographs 
may  be  made  directly.   One  feature  of  this  plotting  routine  to  be  kept  in 
mind  while  examining  the  results,  is  that  it  selects  the  most  suitable  scale 
for  each  specific  plot  within  the  sweep  of  the  beam.   Thus,  if  the  final 
solution  is  calculated  using  different  schemes,  and  if  in  any  one  solution 
so  calculated,  there  is  a  peak  or  overshoot,  the  plotting  subroutine  auto- 
matically changes  the  scale  for  this  particular  curve.   If  now  the  different 
final  solutions  were  photographed  on  a  single  film,  the  solution  with  the 
overshoot  would  be  displaced  from  the  other  solutions,  even  though  the  curves 
represented  the  same  numerical  values,  except  for  the  overshoot. 
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6.   RESULTS  AND  DISCUSSION 

6.1  Comparison  of  the  Different  Schemes 

Figures  k{&)   and  Mb)  show  the  computed  solution  with  boundary- 
conditions  I  and  II,  respectively.   Figure  h(a)   contains  all  the  five  schemes 
tested.   The  four  explicit  schemes  are  practically  undistinguishable  as • they 
overlap  over  one  another  in  the  photograph.   The  smoother,  less  accurate 
curve  is  the  solution  from  the  Crank-Nicolson  scheme.  Figure  k(h)   contains 
only  the  solution  from  the  four  explicit  schemes.   The  mean  square  difference 
of  the  solution  computed  with  the  explicit  schemes  was  only  0.008  compared 
to  an  ideal  step  function,  while  for  the  implicit  scheme  it  was  0.04.   As  can 
he  seen  from  the  figure,  the  accuracy  of  all  four  explicit  schemes  compares 
faborably.   The  implicit  scheme  contains  an  additional  iteration  for  each 

computed  time  step  in  the  evolution  of  the  steady  state  solution.   The 
overall  accuracy  of  the  implicit  scheme  is  dependent  on  the  accuracy 
of  the  intermediate  iterative  step  as  well  as  the  accuracy  of  the 
implicit  scheme  itself.   The  intermediate  iterative  step  also  consumes  a 
good  deal  of  computational  time,  making  the  implicit  scheme  slower  than  the 
explicit  schemes. 

The  Brailovskaya  scheme  yields  results  comparable  in  accuracy  to 
the  other  explicit  schemes  provided  the  time  step  is  well  within  the  stability 
limits  dictated  by  the  linear  analysis.   This  limitation  on  the  time  step  is 
even  more  stringent  with  the  Lax-Wendroff  scheme.   The  Cheng-Allen  scheme, 
developed  from  the  Brailovskaya  scheme,  possesses  the  best  stability  properties. 
Its  theoretical  stability  limit  is  independent  of  the  Reynolds  number  based 
on  the  linearized  analysis.   With  this  scheme  a  larger  time  step 
can  be  used  at  low  Reynolds  numbers  than  is  permissible  with  the  other 
explicit  schemes.   It  was  found  by  experiment  that  even  at  high  Reynolds 
numbers,  where  the  C.F.L.  criterion  is  less  than  the  viscous  limit,  a  much 
larger  time  step  could  be  used  with  the  Cheng-Allen  scheme. 
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Figure  4(a) 


H 1 1 1 1 h 


FIGURE  4(b) 


FIGURE  4.    ALL    FIVE    SCHEMES 
RE  =60 

dt=dxLre_ 

4 
DX=  1/20 
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Tests  with  the  Lax-Wendroff  scheme  show  that,  as  the  high  Reynolds 
number  inviscid  limit  is  approached,  this  scheme  yields  the  most  accurate 
results.   Figure  5  shows  the  computed  solution  at  a  Reynolds  number  of  1000. 
Two  curves  are  distinguishable.   The  plot  with  the  maximum  waviness  is  the 
solution  obtained  from  the  Brailovskaya,  Cheng-Allen  and  Dufort-Frankel 
schemes  superimposed  on  one  another.   Within  this  curve  is  another  curve, 
obtained  from  the  Lax-Wendroff  scheme  which  is  smooth  except  for  wiggles  at 
either  end  of  the  discontinuity  at  x  =  0.5*  All  other  conditions  being -the 
same,  the  Lax-Wendroff  scheme  performs  better  at  high  Reynolds  numbers.   A 
stable  formulation  of  the  Lax-Wendroff  scheme,  for  viscous  flows  is  however 
quite  difficult.   Before  the  stable  scheme  (^-2),  (^3)  "was  arrived  at,  three 
other  formulations  were  experimented  with  but  without  success.   A  similar 
difficulty  has  been  experienced  by  Rubin  and  Burstein  [5]. 

Table  I  shows  the  number  of  iterations  and  the  computational  time 
(sec)  for  a  typical  run  with  Reynolds  number  Re  =  60,  and  the  number  of  space 
intervals  =  k0. 


2 
A  Re    ,  .,     A 
At  is  equal  to  — r —  and  At,-  =  Ax, 


with  three  intermediate  steps  in  between.   Schemes  1  through  5  correspond  to 

the  Brailovskaya,  Dufort-Frankel,  Cheng-Allen,  Crank-Nicolson  and  Lax-Wendroff 

schemes.   The  number  of  iterations  for  the  Crank-Nicolson  scheme  is  not 

given,  because  the  Crank-Nicolson  scheme  is  implicit  with  an  intermediate  iteration 

step  involved  between  two  successive  time  steps.   This,  in  effect,  is  tantamount 

to  altering  the  scale  of  the  time  step,  and  hence  a  meaningful  comparison 

cannot  be  made. 

The  table  shows  that  the  Dufort-Frankel  scheme  requires  the 
least  number  of  iterations  for  convergence,  and  also  that  the  computational 
time  is  minimum  in  this  case.   This  is  because  the  Dufort-Frankel  scheme  is 
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H 1 


FIGURE  5.    INFLUENCE   OF   HIGH  REYNOLDS   NUMBER 

RE  =  1000     DT  =  001     DX  =  1  /50 

THE   LAX-WENDROFF  SCHEME   IS   THE   BEST  SCHEME 


H 1 1 


FIGURE  6.  INFLUENCE  OF  SPACE   STEP  SIZE 
CHENG -ALLEN  SCHEME 

RE  =  60      DT  =  DX2RE/4 

DX=  1/10,  1/20,  1/40  AND  1/50 
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TABLE  I 


Number  of  Iterations 


Scheme 


At. 


At, 


At, 


At, 


Atr 


129 

116 

13^ 

96 

81 

103 

79 

61 

86 

N.C. 

^7 

75 

N.C. 

37 

68 

121 

88 
N.C. 
N.C. 
N.C. 


Computing  Time    (sec) 


Scheme 


At. 


At, 


At. 


At, 


Atr 


1 

2 

3 

k 

5 

28 

13 

20 

59 

15 

21 

9 

15 

^5 

11 

18 

7 

13 

38 

N.C 

N.C. 

5 

11 

32 

N.C 

N.C. 

k 

10 

26 

N.C 

Re   =   60,      Ax   =   1/1*0 

1  =  BRLA 

2  =  DUFL 

3  =  CH-AL 
k   =   CR-NL 

5   -  LAX-WENDROFF 

N.C.    =  DID  NOT  CONVERGE 
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TABLE    II 

2 

Re   =   60  CHENG-ALLEN  SCHEME  t=        , 


Space   Increment   -   Ax  Computing    time,    sec, 
1/10  1.0 

1/20  2.0 

1/30  8.0 

1/40  20.0 

1/50  29.0 
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a  one  step  scheme  and, therefore,  requires  a  minimum  number  of  arithmetic 
operations  for  the  calculation  of  a  mesh  point.   From  the  table  it  is  also 
seen  that  the  Crank -Nicol son  scheme  requires  the  longest  computational  time. 
At  larger  values  of  the  time  step,  the  Brailovskaya  and  Lax-Wendroff  scheme 
do  not  converge.   They  therefore  require  strict  enforcement  of  the  stability 
criterion. 

6.2  Effect  of  Space  Increment  on  Accuracy  and  Resolution 

Since  truncation  errors  constitute  the  major  source  of  inaccuracies, 
reducing  the  mesh  size  by  reducing  the  space  and  time  increments  should 
result  in  an  improvement  in  accuracy.  Figure  6  shows  the  solution  computed 
with  the  Cheng-Allen  scheme  with  the  space  increment  varying  from  l/lO  to 
1/50. 

With  Ax  =  l/lO  the  solution  is  found  to  be  quite  poor  and  the  final 
solution  inaccurate.   With  Ax  =  l/20  (the  curve  above  that  for  Ax  =  l/l0)s 
there  is  a  substantial  improvement  in  accuracy.   However  with  a  further 
reduction  in  Ax  to  l/30,  l/40  and  l/50,  there  is  no  discernable  improvement. 
Thus  there  is  an  optimum  Ax,  beyond  which,  if  the  space  increment  is  reduced, 
there  is  nothing  to  be  gained  either  in  accuracy  or  resolution. 

A  comparison  of  the  computational  times  for  satisfying  the  same 
convergence  criterion  is  given  in  Table  II.   It  shows  that  for  Ax  =  1/10 
convergence  is  extremely  rapid,  but  the  figure  shows  the  resolution  is  un- 
satisfactory.  For  Ax  =  l/30^  "the  accuracy  is  as  good  as  for  Ax  =  l/50,  but 
the  computing  time  for  reaching  the  steady  state  is  only  one-fifth  of  the 
computational  time  at  Ax  =  1/50. 

Therefore  for  any  given  problem  an  optimum  mesh  size  should  be 
selected  which  gives  the  best  resolution  at  a  maximum  saving  of  computational 
time. 

6.3  The  Influence  of  the  Time  Step 

Tests  with  the  selected  schemes  show  that  the  size  of  the  time  step 
is  governed  mainly  by  stability  considerations.  It  is  found  that  for  a  given 
increment,  the  use  of  a  smaller  time  step  does  not  increase  the  overall 
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accuracy  of  the  steady  state  solution.   See  Figure  7-   The  important  criterion 
is  that  the  size  of  the  time  step  in  obtaining  the  asymptotic  steady  state 
solution  should  be  such  as  to  yield  a  convergent  solution.   Some  schemes  like 
the  Cheng-Allen  scheme  permit  the  use  of  time  steps  in  excess  of  the  time 
increment  demanded  by  stability  criteria.   For  other  scemes,  like  the 
Brailovskaya  and  Lax-Wendroff  schemes,  the  time  increment  should  be  well 
below  the  theoretical  stability  limit  for  convergence. 

In  obtaining  the  asymptotic,  steady  state  solution,  the  cumulative 
effect  of  errors  must  be  considered.  For  example,  if  two  different  time 
increments  At  and  Atp  are  used  such  that  Atp  =  2Atn ,  it  is  found  that  the 
mean  square  difference  of  the  computed  solution  from  the  true  solution  changes 
by  about  0.0001.   Since  the  total  time  T  to  attain  the  steady  state  solution 
is  constant,  the  number  of  iterations  in  time  Np  ~  N  / '  •   The  cumulative 
Fourier  components  of  the  error  function  over  Np  ~  N,  /p  iterations  is  less 
than  over  N-,  iterations.   But  the  amplitudes  are  larger  over  Np  iterations 
than  with  N-,  iterations.   These  two  factors  tend  to  even  out,  producing  the 
same  overall  error  with  Np  iterations  with  large  At  as  with  E,  iterations 
with  small  At. 

In  actual  practice  it  would  therefore  be  advisable  to  select  as 
large  a  time  step  as  permissible  to  obtain  good  resolution  while  ensuring 
stability  of  computation. 

6.k     Boundary  Conditions 

Figures  8  and  9  show  solutions  obtained  with  boundary  conditions 
I  and  II  respectively  under  identical  run  conditions,  and  Table  III  shows  the 
number  of  iterations  and  computational  times  for  reaching  the  steady  state 
also  under  identical,  but  different  run  conditions.  When  the  boundary 
conditions  are  such  that  a  discontinuity  exists  at  one  boundary,  the  number 
of  iterations  and  computational  time  are  both  higher.   This  is  due  to  the 
fact  that  with  B.C.I,  the  effects  of  the  discontinuity  at  the  center  are 
compensated  by  positive  and  negative  effects  on  either  side,  while  with  B.C. 2, 
the  effects  of  the  one-sided  differencing  across  the  discontinuity  take 
longer  to  be  damped  out. 
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FIGURE  7.   INFLUENCE   OF   TIME  STEP.  CHENG-ALLEN  SCHEME 
RE  =  60      DX  =  l/30       DT  =  DX   AND    DT  =  DX2RE/4 
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FIGURE   8.      B.C.I 


FIGURE   9       B.C.IL 

FIGURES  8  AND  9.    INFLUENCE   OF    DIFFERENT    BOUNDARY 
CONDITIONS    ON    EXPLICIT    SCHEMES 
RE  =  60       DX  =  l/30       DT=DX2RE/4 
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Table  III 


Scheme 

1 
2 
3 
5 


Number  of  Iterat 

ions 

B.C.I 

B.C. 2 

65 

96 

52 

76 

63 

90 

58 

86 

Re  =  40 
Dx  =  1/20 
At  =  Ax2Re/4 


Computing  Time 

Scheme 

B.C.I 

B.C. 2 

1 

6 

10 

2 

2 

4 

3 

k 

6 

5 

3 

5 

6.5  Formulation  of  the  Boundary  Conditions 

The  calculation  of  the  first  mesh  point  next  to  the  prescribed 
boundary  requires  some  additional  information.   In  actual  calculations  the 
boundary  points  present  the  greatest  difficulty.   The  results  of  the 
numerical  experiments  with  the  Burgers  equation  stress  the  importance  of 
the  correct  formulation  of  the  boundary  conditions  in  a  way  that  is  as  close 
to  physical  reality  as  possible. 

Tests  were  conducted  with  different  formulations  for  calculating 
the  point  immediately  next  to  the  boundary  with  schemes  1,  2  and  3>  namely 
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the  Brailovskaya,   Dufort-Frankel   and   Cheng -Allen  schemes.      The  Brailovskaya 
and  Cheng-Allen   schemes   need  one   extra  point  beyond  the  boundary   if  the 
formulation  for   the  other   interior   points   and  for   the  point   next  to  the 
boundary   is   to  be  the   same. 


i 

V 

E 

D 

A 

B 

C 

G  1 

*1  2  FIGURE  10 

The  most   common  procedure   in  practice  under  these  conditions    is 
to   assume  reflection  conditions   at  the  boundary,    i.e.    to   stipulate  that   the 
values   of  C  and  A   in   Figure  10   are  the   same.      This   arbitrary   specification 
of  values   leads   to   inaccuracies   and  longer   computational   times. 

Table   IV   compares   the  number   of   iterations   and   computational   times 
for   the  Brailovskaya  and  Cheng-Allen  schemes   for  two  different   formulations 
of  the  same  boundary  conditions  given  by   Set    I.      In  formulation  1  reflection 
conditions   are  assumed  as    shown   in   Figure  10.      Formulation  2   is   arrived  at    in 
a  more  physically  realistic   way. 

Since  the  Burgers   equation   is   analogous   to  the  momentum  conserva- 
tion equation,    the  finite  difference   scheme  for   the  cell  adjacent   to   the 
boundary   is   derived   from  the   integral   formulation  applied  to  this    cell.      Thus 


/"* 


2 


X2    9(u2/2 


dx 


(U8 


Table  IV 

Number  of  Iterations 

Computing  Time 

Scheme 

Formulation  1  Formulation  2 

Formulation  1  Formulation  2 

1 

89            37 

13            h 

3 

UU            38 

h                       3 

Re  =  U0 

Dx  =  1/20 

Dt  =  0.01 
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The  calculations  show  that  convergence  is  more  rapid  and  accurate  with 
formulation  2  than  with  formulation  1.   Hence  in  actual  calculations  the 
boundary  conditions  should  he  formulated  in  such  a  way  as  to  represent  the 
physical  boundary  conditions  as  closely  as  possible  without  introducing 
extraneous  numerical  errors  from  extrapolation  formulae  [19]- 

6.6  Perturbation  of  Boundary  Conditions 

To  study  the  influence  of  arbitrary  perturbations  in  boundary 
conditions,  such  as  would  be  the  case  where  unwanted  and  undetected  errors 
crept  into  an  actual  computation  at  the  boundary,  tests  were  run  where  an 
error  of  fO.OOl  was  introduced  at  the  left  boundary,  and  an  error  of 
-0.001  was  introduced  at  the  right  boundary  with  boundary  condition  I,  at 
every  10th  point  on  the  time  axis.   No  convergence  was  obtained  with  any 
of  the  schemes  when  the  perturbation  was  introduced  at  every  5th  point. 
With  the  boundary  values  perturbed  at  every  tenth  point  convergence 
was  obtained  at  low  values  of  the  time  increment  with  all  the  schemes, 
except  the  Dufort-Frankel  scheme,  Table  V.   As  the  time  increment  was 
increased,  the  Cheng-Allen  and  Crank-Nicolson  schemes  retained  their 
stability  while  the  Brailovskaya  and  Lax-Wendroff  schemes  became 
unstable. 

6.7  The  Convergence  Criterion  and  Convergence  Rates 

Two  convergence  norms  were  monitored  during  the  calculations 
namely 

I  r-n    ni 
max  u .  -  u . 

j     J     J 

and 

J 

v  r  /  nflN2   /  ns2  , 

L  [  (u   )   -  (u  )   ] 

3=1 

Figure  11  shows  how  the  mean  square  difference  between  successive  time 
iterations  varies  as  a  function  of  real  time  plotted  on  the  x  axis.   It  is 
seen  that  the  curve  consists  of  two  sections — an  initial  segment  where  the 
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Table  V 


Boundary  Values  Perturbed  at  Every  Tenth  Point 


Number  of  Iterations 


Scheme 


At- 


At2 
At3 

Ate 


1 

2 

3 

k 

5 

129 

unstable 

136 

127 

97 

tt 

107 

unstable 

unstable 

tt 

88 

tt 

it 

tt 

78 

tt 

tt 

tt 

69 

ti 

Computing  Time  (sec) 


Scheme 

At 

At2 

At3 

At^ 

Ate 


1 

2 

3 

k 

5 

29 

unstable 

20 

59 

16 

22 

tt 

16 

^5 

unstable 

unstable 

tt 

13 

38 

tt 

it 

tt 

12 

32 

ti 

tt 

tt 

10 

26 

tt 
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mean  square  difference  remains  almost  steady  and  is  of  the  order  of  0.1  and 
a  second  segment  where  the  curve  dips  downward  until  it  finally  attains  a 
value  of  order  10"  •   These  two  sections  of  the  curve  depict  two  different 
stages  in  the  evolution  of  a  convergent,  asymptotic,  steady-state  solution 
from  a  given  initial  solution. 

During  the  early  part  of  the  calculation,  the  time  dependent 
solution  is  tending  towards  the  steady  state  solution  at  a  rapid  rate.   In 
fact*  it  is  desirable  to  have  this  segment  of  the  curve  rise  rapidly  if  one 
is  interested  in  the  final  steady  state  solution  alone.   This  would  ensure 
a  fast  rate  of  convergence.   If  the  study  of  the  temporal  development  of 
the  flow  is  the  object  of  the  calculation,  it  would  then  be  preferable  to 
have  this  part  of  the  curve  drop  gradually  all  the  way. 

The  second  segment  of  the  curve  in  Figure  11  shows  a  uniform 
drop  to  a  low  value.   This  portion  of  the  curve  is  significant  from  a 
stability  viewpoint.   The  solution  becomes  unstable  unless  the  mean  square 
difference  tends  to  zero.   Figure  12  shows  the  variation  of  the  mean  square 
difference  in  a  divergent  situation.   It  is  seen  that  although  the  curve 
remains  flat  for  a  major  part  of  the  time,  it  finally  yields  increasing 
values  of  the  mean  square  difference  at  the  long  time  limit,  thus  leading 
to  a  non-converging  solution. 

Figure  13  shows  the  effect  of  the  space  increment  on  convergence 
rates.   The  space  increment  decreases  from  l/lO  through  l/20,  l/40  and  l/50 
as  the  curves  progress  toward  the  right.   With  Ax  =  1/10,  the  initial 
segment  of  the  curve  rises  sharply  and  falls  steeply  to  yield  a  high 
convergence  rate.   A  similar  trend  is  observed  for  Ax  =  l/20,  l/40  and 
l/50  although  the  convergence  rate  is  less  rapid.   The  fast  convergence 
rates  of  the  Cheng-Allen  scheme  render  it  suitable  for  computing  the 
asymptotic  steady  state  solution  without  great  attention  to  temporal 
resolution,  while  the  Dufort-Frankel  scheme  focuses  on  the  temporal  evolu- 
tion of  the  final  solution.   Further  the  Dufort-Frankel  scheme  is  second 
order  accurate  in  time,  while  the  Cheng-Allen  scheme  is  only  first  order 
accurate. 
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FIGURE  11.    DUFORT-FRANKEL    SCHEME 
RE  =  50     DX  =  l/20     DT  =  0.02 


FIGURE  12.    BRAILOVSKAYA    SCHEME 
RE  =30     DX  =  30     DT=0.02 


FIGURE   11  AND  12.    R.M.S.  DIFFERENCE   BETWEEN   SUCCESSIVE 
ITERATIONS    PLOTTED  AGAINST    TIME  TO   SHOW  THE 
CONVERGENCE   CHARACTERISTICS  OF  DIFFERENT  SCHEMES 
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FIGURE  13.    INFLUENCE   OF  SPACE   STEP  INCREMENT  ON  CONVERGENCE 

CHARACTERISTICS.    CHENG -ALLEN    SCHEME 

RE  =  60      DT=DX2RE/4     DX=1/10,   1/20,  1/40  AND  1/50 
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Figure  li+  shows  the  value 

i.e.,  the  difference  between  the  true  and  computed  solutions  plotted  along 
the  y  axis  for  B.C.I.   The  solution  is  quite  good,  except  at  the  discon- 
tinuity where  large  gradients  exist.   This  leads  to  large  values  for  the 
derivatives 


du   d  u   du 

ax'    ax2'    c*3 


at  the  discontinuity.   The  influence  of  this  is  discussed  in  the  next  section, 

6.8.   Influence  of  Reynolds  Number 

Figure  15  shows  the  steady  state  solutions  for  B.C.I  for  Reynolds 
numbers  ranging  from  20  to  200,  calculated  with  the  Dufort-Frankel  scheme. 
Similar  calculations  for  B.C. II  obtained  with  the  Brailovskaya  scheme  are 
shown  in  Figure  l6.   In  either  figure  the  curves  with  the  discontinuity 
spread  out  over  the  largest  number  of  meshes  correspond  to  the  lowest  value 
of  the  Reynolds  number.   As  the  Reynolds  number  increases,  the  discontinuity 
thickness  becomes  smaller  and  smaller.   However,  as  the  discontinuity  becomes 
sharper,  wiggles  form  at  the  transition  as  shown  in  Figure  15 
for  Re  =  200.   This  oscillation  at  the  discontinuity  is  caused  by  the 
phenomenon  mentioned  at  the  end  of  the  last  section.   With  an  increase  in 
Reynolds  number,  the  discontinuity  becomes  stronger,  with  the  result  that 
the  gradients 

cfri   c^  u   cru 
^x'   ^2'   cbc3 


increase  in  magnitude.   Since  we  are  differencing  across  the  discontinuity 
without  making  any  special  provision  for  it,  the  truncation  errors  which  are 
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FIGURE   14.     DIFFERENCE    BETWEEN    TRUE   AND   COMPUTED    STEADY 
STATE    SOLUTION.     BRAILOVSKAYA    SCHEME. 
DX  =  l/20,     DT  =  0.02,      RE=50. 
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FIGURE   15.     DUFORT-FRANKEL    SCHEME 
RE  =  20,   40,  60,  100,  200 
DX  =  l/30      DT  =  DX2RE/4 


FIGURE  16.     BRAILOVSKAYA    SCHEME 
RE  =10,   20,  40,  60 
DX=l/20     DT=DX2RE/4 

FIGURES   15  AND  16    SHOW   THE    INFLUENCE   OF   REYNOLDS   NUMBER 
ON    THE    COMPUTED   SOLUTION. 
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proportional  to  these  derivatives  become  significant  at  the  discontinuity, 
leading  to  the  formation  of  wiggles.   This  is  verified  in  Figure  IT  which  is 
the  same  case  as  Figure  16  hut  with  Ax  =  1/50  instead  of  Ax  =  1/30.   No 
waviness  appears  at  the  discontinuity  at  Re  =  200  because  of  reduced 
truncation  errors,  even  though  we  are  differencing  across  a  discontinuity. 

Conclusions 


Tests  with  the  non-linear  Burgers  equation  have  illustrated  a 
number  of  important  aspects  of  the  finite  difference  solution  of  partial 
differential  equations.   Comparison  was  made  between  explicit  and  implicit 
schemes,  and  also  between  one-step  and  two-step  schemes.   The  tests 
show  that  the  Dufort -Frank el  one -step  scheme  and  the  Cheng-Allen  two-step 
scheme  both  possess  definite  advantages  compared  to  the  other  schemes, 
namely  the  Brailovskaya,  Crank-Nicolson  and  Lax-Wendroff  schemes.   The 
Cheng-Allen  scheme  possesses  the  best  stability  properties  when  the 
boundary  conditions  are  perturbed.   This  allows  the  use  of  a  large  time  step. 
The  Dufort -Frank el  scheme  possesses  good  temporal  resolution,  and}in  addition 9 
is  only  a  one-step  scheme.   It  requires  the  shortest  computing  time  while 
providing  results  which  compare  with  the  Cheng-Allen  scheme  in  accuracy. 

The  size  of  the  space  increment  is  governed  mainly  by  how  good  the 
resolution  should  be.   The  use  of  as  large  a  space  increment  as  possible 
for  the  kind  of  resolution  required  is  recommended.   This  results  in  a 
lower  number  of  mesh  points  and  also  permits  the  use  of  a  larger  time  step. 
Discontinuities  in  the  flow  must  be  properly  treated  and  boundary  conditions 
must  be  as  close  to  physical  reality  as  possible  for  meaningful  results. 
The  use  of  a  large  time  step  results  in  much  saving  in  computational  time. 
For  any  one  problem  an  optimum  balance  among  all  these  factors  must  be 
struck. 
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H 1 1 h 


FIGURE  17.      INTERACTION   OF    SPACE    STEP   SIZE    AND  REYNOLDS 
NUMBER    EFFECTS. 

RUN   CONDITIONS    SAME   AS   FIGURE  15,   EXCEPT    DX=l/50 
INSTEAD  OF   1/30. 
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Application  to  Parallel  Processing  Machines 

The  Cheng-Allen  scheme  requires  data  from  five  points  at  time  level 
n  for  the  calculations  of  one  point  at  time  level  (n+l)  if  the  calculations 
are  made  at  one  pass.   However,  by  calculating  all  the  intermediate  quanti- 
ties at  one  pass,  and  then  computing  the  required  final  values  at  time 
level  (n+l)  at  the  second  pass  through  the  arithmetic  units,  the  calcula- 
tion of  a  point  at  level  (n+l)  would  always  require  only  three  values  stored 
at  different  locations  at  time  level  n.   The  Dufort-Frankel  scheme  requires 
fewer  arithmetic  operations,  resulting  in  a  saving  in  computational  time 
especially  in  three  dimensional  problems.   Should  storage  space  "be  at  a 
premium,  the  values  at  time  level  (n+l)  can  "be  calculated  at  one  pass  from 
stored  values.   Thus  both  the  Cheng-Allen  and  Dufort-Frankel  schemes  possess 
desirable  characteristics. 
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